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ABSTRACT 

We discuss a new algorithm to generate multi-scale initial conditions with multiple levels 
of refinements for cosmological "zoom-in" simulations. The method uses an adaptive con- 
volution of Gaussian white noise with a real space transfer function kernel together with 
an adaptive multi-grid Poisson solver to generate displacements and velocities following 
first (ILPT) or second order Lagrangian perturbation theory (2LPT). The new algorithm 
achieves RMS relative errors of order 10~^ for displacements and velocities in the refine- 
ment region and thus improves in terms of errors by about two orders of magnitude over 
previous approaches. In addition, errors are localized at coarse-fine boundaries and do not 
suffer from Fourier-space induced interference ringing. An optional hybrid multi-grid and 
Fast Fourier Transform (FFT) based scheme is introduced which has identical Fourier 
space behaviour as traditional approaches. Using a suite of re-simulations of a galaxy 
cluster halo our real space based approach is found to reproduce correlation functions, 
density profiles, key halo properties and subhalo abundances with per cent level accuracy. 
Finally, we generalize our approach for two-component baryon and dark-matter simula- 
tions and demonstrate that the power spectrum evolution is in excellent agreement with 
linear perturbation theory. For initial baryon density fields, it is suggested to use the local 
Lagrangian approximation in order to generate a density field for mesh based codes that 
is consistent with Lagrangian perturbation theory instead of the current practice of using 
the Eulerian linearly scaled densities. 



Key words: cosmology: theory, large-scale structure of Universe 
methods: numerical 



galaxies: formation 



o 
o 



1 INTRODUCTION 

Cold dark matter density perturbations drive all of cosmologi- 
cal structure formation (e.g. |Peebles|1982l [Davis et al.|1985) in 
the current paradigm. Their non-linear growth is followed with 
numerical simulations over an enormous range of scales. From 
scales of hundreds of Megaparsecs (e.g. jSpringel et al.|[2"005| 
[Pichon et al. 2010) studying cosmological large-scale struc- 



ture; to several hundred kilo parsecs modeling the formation 
and evolution of galaxies including aspects of the interstellar 
medium, molecular cloud-formation and black-hole accretion 
processes (e.g.|Di Matteo et al.|20Q8l|Wise Abel|2008||Cev 



erino & Klypin 2009|| 'Bournaud et al. 2010); and even the 
formation of the first proto-stars formed at scales of a few so- 
lar ra dii (e.g. [Abel et al.|2002|[Yoshida et al.|2008| [Turk et al 
2009) are studied in cosmic context. 



The fluctuations in the standard model would extend from 
solar system size to giga parsec scales. No single calculation 
would be able to capture all these scales. One may therefore 
strive to carry out simulations in as large volumes as possible 
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to study a range of representative environments for the for- 
mation of particular cosmological objects. At the same time, 
however, high resolution is necessary to sample the underlying 
gravitational potential wells. A popular way to achieve large 
simulation boxes with as high as possible resolution for indi- 
vidual objects is the so-called "zoom-in" technique (following 

Navarro White||1994| . In 



versions of e.g. Katz et al. 11 1994 



this technique, a small region encapsulating the formation his- 
tory of the object of interest is studied with much higher reso- 
lution inside of a lower resolution representation of its larger- 
scale cosmic environment interacting with the object through 
long-range gravitational forces. The aim is to represent both 
the variance of the fluctuations on scales of the object (sample 
variance) due to the use of a large volume which provides a fair 
sample of cosmic environments (cf. eg. Crain et al.|2009' Hahn 
et al. 2010), as well as the smaller scale fluctuations affecting 
the structure of the object directly. 

Typically, two approaches have been followed for initial 
conditions of "zoom-in simulations" . Either, initial conditions 
are generated on a uniform grid at the full resolution (which 
can amount to a huge need for computational resources for 
deeply nested grids, see e.g. [Prunet et al.||2008l [Stadel et al.| 
2009) and subsequently degraded by either averaging over 
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regions of 2^^ particles, where d is the de-refinement factor 
outside the high resolution region, or by Fourier resampling. 
There are several disadvantages with this approach: First, the 
computational requirements are immense compared to the so- 
lution sought; and second, the simple averaged velocity fields 
are neither continuous nor differentiable across coarse-fine 
boundaries, leading to spurious shocks in baryonic simulations. 
Alternatively, e.g. the Grafic-2 code ( Bertschinger 2001) is 
used which employs discrete Fourier transforms to add small 
scales perturbations to coarser perturbations. The Grafic- 
2 approach however requires the use of an anti-aliasing filter 
which damps out small-scale power and shows oscillatory er- 
rors at the few per-cent level that result from the combination 
of discrete Fourier transforms of different resolution. Further- 
more, it does not incorporate the sampling of the real-space 



transfer function advocated by [Pen (1997) and Sirko (2005). 
This latter aspect is a feature of all codes that generate initial 
conditions by an inverse Fourier transform of a sampled power 
spectrum on a three-dimensional lattice and thus also plagues 



most applications of the previous approach. Sirko (2005) has 



shown conclusively that this conventional method can lead to 
a significant error in the real-space statistical properties, in 
particular for small cosmological box sizes. 

This paper extends on previous work by proposing a new 
algorithm to generate multi-scale nested initial conditions (as 
[Bertschinger^ ^200 1 ) , but using a real-space convolution ker- 
nel to represent the Fourier space transfer function of density 
perturbations (cf. [Salmon] [19961 |Pen||1997| |Sirko||2005t in a 
multi-scale convolution algorithm combined with an adaptive 
multi-grid algorithm to determine initial velocities, particle 
positions and densities compatible with Lagrangian cosmolog- 
ical perturbation theory up to second order. More specifically, 
the approach uses mass conservation constraints to determine 
the convolution kernel from an accurate real space represen- 
tation of the transfer function kernel on a hierarchy of levels, 
which are used in FFT convolutions to obtain the density field 
from a hierarchical white noise field, and a two-way interface 
multigrid Poisson solver that achieves smooth gradients across 
coarse-fine boundaries. The resulting errors are confined to the 
boundaries and at the level of a few 10~^ in units of the stan- 
dard deviation of the respective fields for the interior of the 
refined region. 

We furthermore show that using a hybrid Poisson solver 
(using multi-grid for the hierarchy of grids and FFT for the 
finest grid), our method produces particle distributions that 
closely match the input power spectrum also on the small- 
est scales and make our approach attractive also for uni-grid 
simulations. All the methods discussed in this paper are im- 
plemented in the initial conditions code MusiC (MUlti-Scale 
Initial Conditions). The code also includes an implementa- 
tion of a linearized Boltzmann solver based on Linger fMa" 



& Bertschinger 1995) that allows the calculation of cosmolog- 



ical density and velocity perturbations for dark matter and 
baryons. 

The paper is organised as follows. In Section [2] we dis- 
cuss a multi-scale convolution method that produces a nested 
density field that is consistent with both the input correlation 
function and the input power spectrum, and which is then used 
as the source field for Lagrangian perturbation theory. Com- 
puting displacement and velocity fields in Lagrangian pertur- 
bation theory (LPT) requires a numerical solution to Poisson's 
equation. In Section [S] after briefly summarizing LPT, we de- 
scribe the adaptive multi-grid algorithm that we use to obtain 



displacement and velocity fields. We also introduce a hybrid 
Fourier-space/multi-grid based Poisson solver which has bet- 
ter Fourier-space properties on small scales compared to the 
pure multi-grid algorithm. In Section |4] we discuss the perfor- 
mance of our method compared to initial conditions generated 
purely in Fourier-space both in unigrid and nested grid set- 
ups. We generalize our approach to initial conditions for two- 
component dark-matter and baryon simulations in Section [5] 
and propose the usage of the local Lagrangian approximation 
to generate inital baryon density fields for grid codes before 
we summarize and conclude our results in Section [G] 



2 GENERATING DENSITY PERTURBATIONS 

In this section, we summarize methods to generate Gaussian 
random fields that follow a prescribed power spectrum and 
act as source terms for density and velocity perturbations in 
Lagrangian perturbation theory. We give particular attention 
to a convolution kernel based approach which has favourable 
properties in a multi-scale nested grid set-up. 



2.1 Computing the seed density field 

Consider an over-density field ^(r) that is completely described 
by its power spectrum P{k) =< ^(k) ^*(k) >, where the tilde 
denotes the transformed function of a Fourier transform pair. 
The amplitude of density fluctuations is usually expressed in 
terms of the transfer function T(A;), which is defined such that 



(1) 



where Ug is the (constant) power spectrum spectral index after 
inflation, and a is a normalization constant. The first step 
in setting up initial conditions for cosmological simulations 
thus consists in generating a white noise sample of random 
values /i(r) (typically sampled from a Gaussian distribution 
with zero mean and unit variance) and to require that their 
amplitudes follow a specific power spectrum P{k). This can be 
achieved by multiplying the Fourier transformed white noise 
field Ji(k) with the square root of the power spectrum, i.e. for 
all k representable on a grid of given resolution set 



5(k) = VpOWD M(k) = a ||k||"»/^ r(||k||) M(k). (2) 

The real space over-density field ^(r) is then obtained by in- 
verse Fourier transformation, and we call this procedure "/c- 
space sampling". 

It is interesting to note that a product in Fourier space 
simply corresponds to a convolution in real space (cf. also EF] 
stathiou et al.|p^85 Salmon|[l996 ), i.e. equation ([2| is equiv- 
alent to 



5(r)=T(||r||) *Mr), 



(3) 



where T(r) is the real-space counterpart of T{k) = 
ak^''^'^ T{k), and 'V denotes a convolution. 

It is thus mathematically equivalent whether equation Q 
is evaluated in Fourier space, followed by an inverse trans- 
form, or whether equation (|3| is evaluated using an inverse 
transform of T{k) followed by the convolution. Most cosmo- 
logical initial conditions generators follow the first approach 
(e.g |Bertschinger|200l| ), while e.g. [Salmon] ( |1996| ) , |Pen| ( |1997 ) 
and Sirko ( 2005 ) use the second or variations thereof. The dis- 



crete realizations of the density fields derived with the two 
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Figure 1. The real space transfer function T(r) at 2; = 50 for total 
matter density perturbations obtained using the FFTlog method 
from a tabulated k-space transfer function T{k). The transfer func- 
tion is positive where the line is solid and negative where it is dotted. 
The inset shows the baryon oscillation peak in linear scale where 
the change of sign occurs. 



approaches will have significant differences. In particular, em- 
ploying eq. ([2]) forces periodicity of the real space transfer 
function on box scales that leads to an underestimation of the 



two-point correlation function (see Section 2.2 and e.g. Pen 
and Sirko|2005 for a detailed discussion) . This is particu 



1997 



larly relevant for small cosmological boxes (L < 100 Mpc), 
where the periodic contribution is non-negligible. 

The convolution kernel T(r) can be thought of as a "prop- 
agator" of a single Gaussian white noise fluctuation /iq(x) = 
fc(x — q) at location q, where is the Dirac ^-function. The 
time evolution of each single such perturbation is given by 
the time evolution of T(r). Instead of a variation of the tradi- 
tional method of sampling Fourier modes according to a given 
power spectrum, this real-space picture has thus an intuitive 
physical motivation: The convolution operation imprints the 
density perturbations due to the white noise /i(x) onto space. 
Apart from the better behaviour of the two-point correlation 
function, we follow the convolution approach since it allows 
for an easier treatment of nested grids that are separated by 
sharp boundaries in real space. To achieve this, a method will 
be developed that correctly overlaps the "propagator" with 
such coarse- fine boundaries. We give an account of how we 
numerically compute eq. ([3| in the next section and of the 
multi-scale convolution approach for nested grids in Section 
01 



2.2 The real-space transfer function 

Since the transfer function T{k) is typically computed numer 
ically using a spectral Boltzmann solver - e.g. the Boltzmann 
solver included wi th our code, Linger jMa &: Bertschinger 
' 1995| ), CMBfast dSeljak k Zaldarriaga|| 19961 or Cam^ or 



^ Camb, written by A. Lewis and A. Challinor, can be obtained 
from http://camb.info 



taken from a fitting formula (e.g. Eisenstein & Hu 1998) - 
its Fourier transform is not available a priori. Thus, in the 
approach we follow, an accurate real-space transfer function 
Tn{r) is first calculated which is subsequently applied as a con- 
volution kernel to the Gaussian white noise field using FFT 
convolution. Assuming a spherically symmetric transfer func- 
tion T{k), one has 



TTTTa [ ^e^) exp(ix- k)d^fc 
_4n_ /-^ siii(fer) 



dfc. 



(4) 
(5) 



This real space transfer function, Tn{r), can be computed 
using the FFTlog algorithm ( T almaiilfTsiVsl |Hamilton||2000 ) 
from an input k-space transfer function T{k). FFTLOG uses 
the fact that the Fourier transform of a log-log function can 
be written as a Hankel transform which can be carried out by 
FFT. The resulting transform is highly accurate over many 
orders of magnitude. We kindly refer the reader to |Hamilton| 
(2000) for details. The r = component needs to be computed 



separately. It is obtained by a three dimensional numerical 
integration over the range of Fourier modes included in the 
simulation. This r = component sets the white noise level of 
the density perturbations. 

Figure [1] shows the result of the transform from tabulated 
k-space transfer functions at z = 50 for matter density per- 
turbations generated with a ID linear Boltzmann solver (cf. 
Ma k Bertschinger||1995 ). 

Finally, the convolution in eq. (|3| has to be com- 
puted numerically. To this end, we compute a discretisation 
T{yiijk) = 77^(|||xi^■fc|||), with x^^-fc = {ihjh,kh), {i,j,k) G 
[—N/2 + 1, over the entire simulation volume of length 

L = hN. Next, a white noise random field ii(yiijk) is created 
on the same discretisation. Then both T{xijk) and ii{yiijk) 
are transformed by FFT, multiplied element-by-element with 
one another and the result is inverse transformed to yield 
^(xijfc). Note that when sampling on a finite grid, the real- 
space sampled transfer function is no longer spherically sym- 
metric when transformed by FFT to three dimensional Fourier 
space. However, equivalent ly, the traditional /c-space sampled 
transfer function is not spherically symmetric in real space. 

The analysis performed in this paper uses a real-space 
transfer function periodically replicated on the box length. 
This abandons the favourable properties of a truncated cor- 



relation kernel discussed in detail by Sirko (2005|. However, 
this choice allows comparison with the classical approach that 
samples the k-space transfer function. For actual applications 
it is recommended that the truncated transfer functions should 
be used (which is a simple parameter in our code). 

In Figure [2] we show the auto-correlation function of the 
transfer function kernels (equivalent to the two-point corre- 
lation function). As demonstrated by Sirko (2005'), the non- 



periodic real space kernel perfectly reproduces the input two- 
point correlation function. The periodic correlator, for which 
we have not subtracted the mean in this case, underestimates 
the correlation function at large radii and the /c-space sampled 
kernel shows an even stronger drop at large radii. Note that 
the big drop in the /c-space sampled correlation function is 
somewhat spurious. It is partly from the periodic component, 
as the periodic real space kernel shows a similar suppression. 
However, it is mainly due to an additional integral constraint. 
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Figure 2. Auto-correlation function of the transfer function (equiv- 
alent to the two-point correlation function) for /c-space sampled 
(dotted red) and real space transfer functions that are non-periodic 
(dashed green) and periodic on the box length (dot-dash blue). The 
gray line indicates the two-point correlation function determined 
from the input power spectrum at z = 45. All auto-correlation 
functions have been computed for a 500 /i"-*^ Mpc box with 256^ 
resolution. 



For the /c-space sampled kernel, P{k = 0) = 0, so that 



: 



L 



e(x)d^ 



; 4:71 



J ^(r)r^dr, 



(6) 



where V is the simulation volume and the latter equality holds 
in the case of a spherically symmetric correlation function and 
amounts to the usual integral constraint on the two-point cor- 
relation function. This integral constraint however holds for a 
finite T> in the traditional /c-space sampling approach (while 
D ^ in our case) and the correlation function is thus offset 
by an additive constant which is equal to the integral over the 
correlation function outside the simulation volume. This ad- 
ditive constant leads to the additional deviation seen between 
the /c-space sampled correlation function and the periodic real- 
space equivalent which is simply amplified by the multiplica- 
tion with r^. Any box with zero mean overdensity will there- 
fore fulfill the integral constraint over the box instead of over 
an infinite volume leading to a similar discrepancy. Boxes with 
a non-zero mean overdensity would circumvent this problem 
and the difference in mean overdensity can be incorporated 
in a change to different effective cosmological parameters. We 
will not investigate further the possibility of simulations with 
non- vanishing mean overdensity but kindly refer the reader to 



Sirko ( 2005 ) for a detailed discussion of this possibility. 



2.3 Generating a nested initial density field 

In order to generate an accurate multi-scale representation of 
an initial density field, we will now describe our algorithm 
to perform the convolution in eq. ^ on nested grids. Par- 
ticular care is taken to achieve mass conservation, as viola- 
tions would affect the entire domain when the density field is 
used to source the displacement and velocity fields (cf. Sec- 
tion [sj. As in the Grafic2 approach( Bertschinger|2001 ) , this 
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Figure 3. The volume subdivision based nested grid structure used 
in the multi-scale algorithms. Shown are two additional refinement 
levels in which a parent cell can be refined into eight children cells. 
Note that centers of children do not coincide with the parent grid 
cell centres. 



amounts to linear constraints between the levels. However, un- 
like Grafic2, these constraints are not only top-down (i.e. 
the coarse level provides constraints for the fine) , but partially 
bottom-up to achieve a more conservative algorithm and high- 
resolution sampling of the real-space convolution kernel. 



2.3.1 White noise refinement 

In order to generate initial conditions for a nested subdomain 
on the volume subdivision grid structure used in our approach 
(Figure [3|, the white noise for the subgrid has to be consistent 
with that of the coarse grid. Two separate approaches can be 
followed that are locally mass-conserving (in the sense that 
the children cells average to the parent cell locally rather than 
that only the total mass in the subgrid domain is conserved): 

The first approach uses the Hofmann-Ribak algorithm 
( Hoffman &: Ribak||1991 ). First an unconstrained white noise 
field ct;^^^ is generated for level i -\- 1 which has 8 times the 
variance of the noise field for level i in the case of a refine- 
ment factor of 2. Next, for each group of eight children cells 
on the fine grid, the mean is replaced by the respective value 
of the noise field on the coarser level i (cf. also Pen 1997 
|Bertschinger||2001 t. 

Since this approach does not retain the Fourier modes 
present on the coarse grid, another approach can be devised 
that achieves this: (1) again, an unconstrained white noise 
field cj^^"^ is generated for level + 1 at 8 times the coarse 
level variance. (2) In order to preserve Fourier modes that are 
present in the coarse grid an FFT of both the fine grid 
and that part of the coarse grid that overlaps the fine grid is 
performed. Then all modes k ^ '^Ny,i up to the Nyquist wave 
number of the coarse grid are replaced with the respective 
Fourier coefficients from the coarse grid, followed by an inverse 
transform. (3) In order to ensure conservation of mass, the 
Hoffman-Ribak algorithm is applied in reverse, i.e. the coarse 
grid white noise values are replaced by the average over the 
eight children of a coarse cell inside the refinement region. 
Since the Fourier interpolation conserves the mass of the entire 
subgrid, as does the averaging over children, the global mass 
in the simulation domain is also conserved. 

Further levels can be computed by repeating first step (2) 
for all levels at increasing resolution, followed by the correction 
step that replaces all coarse cells at all coarser levels with the 
average over the eight fine cells, starting at the finest level. 
Thus, the Hoffman-Ribak constraint is fulfilled while at the 
same time Fourier modes are preserved between grids. 
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In the following, we will describe how the convolution with 
the real-space transfer function, eq. ([s]), is performed when a 
refinement region is present. 

2.3.2 Generating convolution kernels 

In Figure [4] we show schematically the set-up for one addi- 
tional refinement level. The top grid domain Q consists of the 
domain Q2 covered by a refinement grid and the non-refined 
part Q\Q2- The refinement region Q2 at twice the resolution 
is given by Q^ ^2,p represents the padding around Q2 to twice 
the length per dimension that is needed to perform isolated 
Fourier convolutions. 

The next step in our method is to generate convolution 
kernels T(r) for all required levels. In contrast to the Grafic2- 
approach, we construct these kernels purely in real space, 
starting from the finest level. For the finest level this consists 
of a simple evaluation of the real-space transfer function Tn{r) 
on Q'uQp, i.e. T^(x^^-fe) = 7^(|| |xi^-fc || |) for all Xijk eQ'uQp 
with the origin x = placed at the centre. 

As for the white noise, care has to be taken to main- 
tain local conservation of mass. Thus, rather than sampling at 
coarser resolution, for the next coarser level i — 1, the already 
discretized part of the convolution kernel for the fine-level grid 
needs to be restricted to the coarser grid in a conservative way. 
This can be achieved by averaging the contributions due to all 
eight children cells when convolved with the kernel, i.e. we 
compute 

T'-\x,y,z) = ^ J2 T\x + i,y + j,z + k), (7) 

for those grid points {x,y,z) that correspond to the refined 
part of the grid. This provides the convolution kernel for a 
volume of the size of the refined region including the padding, 
i.e. for Q2 U ^2,p- The remaining cells outside that region, i.e. 
on Qi\(Q2 U ^2,p) are again sampled directly from the real- 
space transfer function: T^~^(x^jfc) = Tn{r) for all x^^fc G 
^i\(^2 U ^2,^)- This last step is in violation with strict mass 
conservation since eq. ([7| does not hold outside the refined re- 
gion. Since most of the mass of the convolution kernel is close 
to the centre, this approximation produces however negligible 
errors and avoids the substantial computational overhead of 
an exact evaluation of 77^(r) over the entire domain at the 
finest resolution. The procedure of averaging and newly sam- 
pling uncovered volume is repeated until the coarsest level is 
reached. 

2.3.3 Noise convolution: The first refinement level 

The density field on the top grid level i is determined as in 
the unigrid case by computing 5^ — T^ i< fi^ on Q with periodic 
boundary conditions automatically satisfied by the FFT. For 
the refined region, several contributions will be co-added: 

1. The coarse grid contribution d^t^vse to the refinement re- 
gion is computed by zeroing f/ on Q2 to obtain /j^q^q^, 
computing the convolution 6{ = i<: fifi^Q^ and interpolat- 
ing si from Q2 to Sctarse- Wc usc tri-cubic interpolation for 
this purpose. 

2. On level i -\- 1 isolated boundary conditions apply, so 
has to be padded to twice its size 2N so that an FFT-based 
convolution is still possible. In order to find the density 
^ttif sub-grid alone, on Qp is set to zero to 



a) top grid b) subgrid 

































< 2N ► 



Figure 4. Schematic representation of a set-up with one refinement 
level. The left panel a) shows the top grid Q which consists of a 
region Q2 covered by a refinement grid and the non-refined part 
r2\r22. The right panel b) shows the sub-grid domain which is 
equivalent to Q2 in the left panel but has twice the resolution. To 
this sub-grid a padding region flp will be added when performing 
the FFT convolution with isolated boundaries and is denoted by 
Cl2,p on the top grid. 

obtain /x^t^, followed by computing the convolution 6^^^^ = 

3. In order to reduce errors at the boundary, a correction 
term 6^'^^ is added that accounts for the fiuctuations just 
outside Q^ To this end, the coarse grid value of the white 
noise field is subtracted from each of the 8 children cells, 
equivalent to "unapplying" the Hoffman-Ribak method. We 
then zero the result of this operation on Q', so that it 
is non-zero only on Qp and obtain Since boundary 
conditions are isolated, it would be necessary to truncate 
the transfer function in order to have a non-periodic un- 
constrained white noise field. We however found that a 
truncation introduces larger errors than assuming period- 
icity (on scales larger than the subgrid). Finally, the FFT- 
convolution S^'^^ = j^^+i ^ is evaluated. 

4. Finally, the result of the previous step is restricted also to 
the coarse grid, i.e. onto ^2,^, in order to include informa- 
tion about fiuctuations on smaller scales. 

Finally, all three contributions are added to find the refined 
density field 

^ - ^self + ^coarse + <>bnd l^J 

on . In order to further reduce errors due to the boundary, 
we allow for optional additional padding of the subregion with 
a few grid cells and cut out the desired region at the end. 

2.3.4 Noise convolution: Further refinement levels 

For further refinement levels, the same steps 1-4 from above 
are repeated at level i + i, i ^ 2, with the only difference 
being that the coarse grid contribution Sctarse is computed us- 
ing isolated boundary conditions on ^ + i — 1, i.e. with a zero 
padded random field Furthermore, all coarse contri- 

butions have to be interpolated down to this level. Let V be 
the interpolation operator, then 

5^+^ = Cf+Cd+ E T^'-'l^^^e], (9) 

j=0,...,2-l 

where [•] indicates a successive interpolation between q lev- 
els. We use a conservative tri-cubic interpolation for this pur- 
pose, i.e. we first perform a normal tri-cubic interpolation and 
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then correct the eight children cells with an additive constant 
to ensure that their average equals the coarse cell value. 



3 INITIAL PARTICLE POSITIONS AND 
VELOCITY FIELDS 

In this Section, we briefly summarize the application of first 
and second order Lagrangian perturbation theory to obtain 
the initial displacement and velocity fields which are based 
on solutions of Poisson's equation. For details, we kindly refer 
the reader to the wide field of existing literature on Lagrangian 
Buchert et al.||1994| [Bouchet et al 



3.1.2 Second order perturbations 



perturbation theory (e.g 
19951 



Scoccimarro 



1998 



Bernardeau et al. 2002). We then 



summarize the multi-grid algorithm which solves Poisson's 
equation numerically before we discuss its extension to the 
adaptive multi-grid algorithm which provides solutions to 
Poisson's equation on nested grids. Finally, we examine several 
methods to obtain velocity and displacement fields that have 
the same behaviour at large wave numbers as those obtained 
with the traditional /c-space sampling. 

3.1 Lagrangian perturbation theory 

3.1.1 First order perturbations 

Lagrangian perturbation theory describes the evolution of 
density perturbations in the rest- frame of the fluid. The posi- 
tion X of a fluid element at time t with respect to its initial 
position q, and the respective fluid velocity, can then be writ- 
ten as 



x(t) = q + L(q,t), x(t) = -L(q,t) 



(10) 



where L(q) we call the "displacement field" which is derived 
using perturbation theory. 

It can be easily shown that at first order in the pertur- 
bations (cf. Zel'Dovich 1970), the displacement field L can be 
written as the gradient of a potential $ which is proportional 
to the gravitational potential 0, 

^^''^ = ~ 3HiJD+{t) - ^+'^^^ ^^^^ 

where Ho is the Hubble constant, a is the expansion factor at 
time t, Dj^(t) is the growth factor of linear density perturba- 
tions (i.e. the time evolution of the growing mode) and (j) is 
the gravitational potential obeying Poisson's equation 



^c<^>{^,t) = ^HW5{^,t). 



(12) 



Since the velocities are given by the gradient of a potential, 
velocities are irrotational, i.e. V x x(t) = 0, in this approxi- 
mation. 

Note that the Gaussian over-density field 5 is the source 
field of the displacements. It is not the density field that is 
measured after displacing the fluid elements, which is non- 
Gaussian. We give a derivation of the latter in Section |5.3| 
The Gaussian field 5 should not be used to impose an initial 
density field for the baryonic component. 

In order to obtain the displacement vectors from the ini- 
tial over-density field ^, Poisson's equation has to be solved 
numerically. The most common approaches use an FFT based 
Poisson solver (e.g. Bertschinger 2001 ), while we chose a multi- 
grid based Poisson solver as it can be easily extended to an 
adaptive multi-grid solver which is able to handle nested adap- 
tive grids in a very natural way. 



cimarro 1998 



y. Munshi et al.|1994 


|Scoc- 


Tatekawa & Mizuno 


(2007) 



that first order Lagrangian perturbation theory (cf. Section 
3.1.1) might not be accurate enough for current simulations 



as it leads to e.g. significantly underestimated higher order 
moments of the density probability distribution functions at 



early times. In particular, Munshi et al. (1994) have shown 
that 2LPT matches the skewness of the density field for top 
hat collapse, while 3LPT would in addition also match the 
kurtosis, and so on. At second order, the displacement field 
contains not only contributions from the gravitational poten- 
tial, but also from a second-order potential that is sourced by 
the off-trace components of the deformation tensor, i.e. 

L(q,t) = D+(t) Vq$(q,t) + D2(t) Vq*(q,t), (13) 

where * obeys the Poisson equation Aq*I'(q, t) = T(q, t), with 

^(q-*) = -\Y. [i^o^^'^f - id,,d,,'i) {d,,d,,^)\ , (14) 

and D-^(t) is the growth factor of linear perturbations, and 
D2{t) ^ |i^+(t). Adding second order perturbation theory is 
thus algorithmically identical to repeating the steps for first 
order displacements: After computing the source- field r using 
finite differences, Poisson's equation can be solved numerically 
using the (adaptive) multi-grid scheme. A similar adaptive 
approach to generate initial conditions at second order has 
been made by Jenkins (2010) who use a tree-PM method to 



evaluate the second order contribution. Note that while our 
method allows for 2LPT, we only use ILPT in this paper to 
aid the comparison with the typical practice in the field. 



3.2 Multi-grid solution of Poisson's equation 

Both first and second order Lagrangian perturbation theory 
for velocity and displacement fields require the numerical so- 
lution of Poisson's equation followed by calculating gradients. 
This can be achieved by use of the multi-grid algorithm (e.g. 



Fedorenko 1961 Brandt 1973). In order to solve Poisson's 



equation 



A(/)(x) = /(x) on domain Q, 



(15) 



with periodic boundary conditions in our case, we cover Q 
with a hierarchy of grids M°, M^, . . . , of respective grid 
spacing h^,h^ , . . . , h'^ , where /h^^^ = 2, i.e. a refinement 
factor of 2 between multi-grid levels (cf. also Figure [sj. 

Define I^~^ as the restriction and I^^^ as the injection op- 



erator. We use the Full Approximation Scheme (FAS - Brandt] 
1977), see also |Trottenberg et al. (2001), to solve the discrete 
form of equation ( 15 ) on grid I given by 



Vu\x) = f[x) for X G 



(16) 



where L is a finite difference approximation to the Laplacian 
(as in Table 111 and is an approximation to (j) on grid M^. 
The residual r^{x) can then be written as r^ = - L^u^. FAS 
then uses these residuals to correct the solution on level i — 1, 
so that (note that operators do not commute) 



T i-l i-1 _ i I J i-lji-l i 

Li U — L/) r -]- Li Ln U 



is obtained. This is equivalent to solving 



(17) 



(18) 
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Order n Laplacian L Gradient G 



exact: d'^ dx 

2: [1-21] 1 [ _i 1 ] 

4: 12 [ -1 16 -30 16 -1 ] ^ [ 1 -8 8 -1 ] 

6: tIo [ 2 -27 270 -490 270 -27 2 ] ^ [ -1 9 -45 45 -9 1 ] 



exact: —k? —ik 

2: -2 [-cos(/c) + 1] -z sin(fc) 

4: - 1 [cos(2/c) - 16 cos(/c) + 15] - § [- sin(2/c) + 8 sin(/c)] 

6: - ^ [-2 cos(3/c) + 27 cos(2/c) - 270 cos(/c) + 245] - ^ [sin(3/c) - 9 sin(2fc) + 45 sin(/c)] 



Table 1. Finite difference stencils in one dimension for the Laplacian L and gradient operators G up to 6th order (top rows) and their 
respective Fourier transforms L and G (bottom rows). 



on level 



1, where 



Order n 



Flux operator F 



£-1 



jt — ij t t 



(19) 



The grid at level I thus provides additional source terms r to 
the equation at level — 1 in addition to the restricted source 
f^~^ = accounting for the non- commutative nature of 

the operators - eq. (19) is just the commutator [L, I^""*^] of 
the Laplacian and the restriction operator. 

Within each level, corrections are propagated using a 
smoothing (or diffusion) scheme S (u^ , /^), typically a Gauss- 
Seidel sweep. 

The complete FAS multi-grid scheme then consists of the 
following algorithm to solve L^u^ — starting with £ — m: 

1. If ^ = 0, set = 0. 

2. Apply z^i smoothing steps, u{ = S (u{-i^ /^), i = 1 . . . z^i 

3. Calculate the residual, = — llu^ 

4. Restrict the residual, = I^~V^ 

5. Restrict the smoothed solution, u^~^ = 

6. Apply the r-correction, = + L^~^u^~^ 

7. Apply FAS scheme recursively to solve L^~^i^^~^ = /^~^ 

8. Correct the solution = ui^ + (u^~^ — l\~^ui^) 

9. Apply V2 smoothing steps, ul = S (ul_i, /^), i — 



■ ^2 



The first step ensures that the mean of the potential van- 
ishes and that the algorithm converges in the case of periodic 
boundary conditions. For non-periodic boundary conditions, a 
direct solution would need to be computed. Throughout each 
cycle we enforce periodic boundary conditions, whenever u 
is changed. The scheme is to be repeated until the norm of 
the residual, computed after step 9, falls below some desired 
threshold. Since we only call FAS once in step 7, our approach 
uses only V- cycles. 

We found excellent convergence, i.e. a reduction of the 
residual by at least one order of magnitude per iteration, for 
all finite difference approximations of the Laplacian that we 
tested (up to sixth order, cf. Table [T]) using the red-black 
Gauss-Seidel method as the smoothing operation S (u^ , /^) 
with z/i = z/2 = 2 sweeps, simple averaging over the 8 child 
cells of one coarse cell as the restriction operation and 
straight injection of the coarse cell value into the 8 child cells 
as the prolongation operation 

The FAS scheme has been adopted as it operates on the 
solution itself on each level, while the standard multigrid op- 
erates on residuals. This allows us to incorporate conservative 
boundary conditions at coarse-fine boundaries more easily (cf. 
Section 3.3). 



1: 


, [ 


1 ] 




3: 


[ -1 15 


-15 1 ] 




5: 


[ 2 -25 245 


-245 25 


-2 ] 



Table 2. Finite difference flux operators for the Laplacian. Con- 
volved with [ — 1 1 ] , these become the respective Laplacians of order 
(n+1). 



3.3 Adaptive multi-grid 

We will now describe the extension of the FAS multi-grid al- 
gorithm described above to additional nested adaptive grids, 
M^^^^^ . . . M^'^^^, covering non-coextensive subdomains 
i — 1, . . . , /c with Q^*^^ C In our case, the M' are simply 
rectangular grids. Two modifications have to be made: First, 
restriction and prolongation only operate on over- 
lapping regions Q.'^ H Q.'^^^ of the domains. The remainder 
is treated as if it resided at the finest level. Second, 
Poisson's equation on additional sub-grids is solved with the 
coarse grid solution acting as a boundary condition for the 
finer level Poisson equation l/^'^u^^'^ — /^^^. The boundary 
condition is realised by adding ghost zones to the sub-grids 
M' to which boundary values are interpolated. 

For these coarse- fine boundaries, we first use polynomial 
interpolation using only coarse grid information parallel to 
the fine grid surface. Using these intermediate values together 
with values inside the fine grid, another polynomial interpola- 
tion step is performed normal to the fine grid surface (similar 
to [Martin Cartwri ght][l99 6t. The order of the interpolat- 
ing polynomial is chosen identical to the order of the finite 
difference scheme for the Laplacian. This serves as a "guess" 
for the solution on the boundary and is thus a weak Dirichlet 
boundary condition. "Weak" as we will not actually use the 
boundary cell solution as the boundary condition. 

In a second step, we apply additional constraints to the 
interpolated ghost zone values so that the coarse flux matches 
the fine flux across the coarse-fine boundary. This procedure 
can be easily understood by rewriting Poisson's equation in 
the following way: 

A<^(x) = V-V<^(x) = /(x) (20) 

can be integrated over one grid cell (in one dimension for no- 
tational simplicity) to yield 

f^Xi-\-h/2 

dx Fci,{x) dx = rrii, (21) 

/xi-h/2 

where rrii — f f{xi) dx' is the mass contained in cell i, and 



rx 

J Xa 
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Fcj) = ^cc0 is the potential flux. Applying the divergence theo- 
rem, we find that 



F^{xi + h/2) - F^{xi - h/2) = rm 



(22) 



i.e. the mass in cell i generates a flux through the cell surfaces 
given by F^. In order for the adaptive multi-grid scheme to 
be conservative, the flux across common surfaces has to be 
identical (since the mass contained in a coarse cell is taken 
to be identical to the mass inside its eight children). This is 
trivially fulfilled for inner cells but care has to be taken at 
outer coarse- fine boundaries. 

The flux operators F are matched to the order of the 
Laplacian and given in Table [2] They are gradient operators 
of order n — 1 on the cell boundaries. When convolved with the 
first order gradient, [ — 1 1 ] (which corresponds to a deter- 
mination of the net flux) , the respective Laplacian operator of 
order n is recovered. Using the flux operator, the 4 fluxes due 
to the fine grid can be evaluated and then subtracted from 
the coarse flux through the same surface element. This flux 
difference is subtracted from the ghost zone values contribut- 
ing to the fine grid flux so that the sum of the fine fluxes now 
matches the coarse flux. A flux- correction is necessary since 
otherwise the sub-grid induces artificial source terms through 
the ghost zone interpolation and the multi-grid algorithm does 
not converge (in the sense that the residual does not vanish on 
all levels). Note that this flux matching can be thought of as 
constraints on the interpolation scheme. We thus obtain mixed 
boundary conditions: a von Neumann boundary condition due 
to the flux matching, which constrains only one degree of free- 
dom for the boundary cells at a single coarse-fine boundary, 
and a subordinate Dirichlet boundary condition, which con- 
strains all degrees of freedom, but is only used to evaluate 
the fluxes on the fine grid (cf. also Miniati k, Colella||2007 ). 
This results in a two-way interface between coarse and fine 
cells. These boundary conditions ensure that gradients across 
coarse-fine boundaries are smooth. Note that this is not the 
case for the one-way interface multigrid solvers currently em- 
ployed in most cosmological simulation codes due to the need 
for adaptive time-stepping. Using this conservative procedure 
we maintain multi-grid convergence (i.e. the residual reduces 
by at least an order of magnitude per iteration) also with an 
arbitrarily deep hierarchy of adaptive sub-grids. 

Note that our approach to use an adaptive Poisson solver 
is in similar spirit as Salmon (1996) who used a tree to com- 
pute displacements and velocities, or | Jenkinsj ( |201Q| who also 
used a tree to compute the second order term for 2LPT. The 
advantage of the multi-grid method is however that it has a 
well controlled residual to the equation to be solved so that 
errors are easily controllable by setting the convergence cri- 
terium in terms of the residual norm rather than by tuning 
opening angles for the tree. The performance of a tree-based 
Poisson solver is particularly problematic for density fields 
with low constrast. Using a tree has the further disadvantage 
that periodic boundary conditions have to be incorporated 
in a hybrid way (e.g. by using FFT for the top grid). This 
leads to an additional source of errors arising from the long- 
range/short-range split. 



3.4 Fourier analysis of the finite difference 
operators 

Operating in real space requires the use of a finite difference 
approximation to the Laplacian "L" and gradient operators 
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Figure 5. Wavenumber dependent attenuation of perturbation am- 
plitudes due to the finite difference approximations of the Laplacian 
and gradient operators. Shown is the relative attenuation with re- 
spect to the exact solution for 2nd (dotted), 4th (dashed) and 6th or- 
der (dash-dotted) approximations in one dimension. The wavenum- 
ber is in units of the Nyquist wave number k^y. 



"G". In Table [T] we give the standard stencil representations 
for the one dimensional versions of these operators. The three- 
dimensional versions can be obtained by subsequent convolu- 
tion of the one-dimensional operators along all three Cartesian 
coordinate axes. In the bottom half, the Fourier transforms of 
the operators are given together with the exact Fourier trans- 
form of the continuous operators. Since a regularly spaced 
mesh is a Dirac comb, due to symmetry reasons, the Fourier 
transform of these operators takes the form of a cosine series, 
in the case of the Laplacian, and of a sine series, in the case of 
the Gradient. These series are approximations to the respec- 
tive continuous and non-periodic transforms of the operators. 
The sine series representing the finite difference gradients van- 
ish at the Nyquist wave number leading to a suppression of 
small scale modes. Furthermore at low order, the finite dif- 
ference approximation leads to an attenuation also at slightly 
large scales. We will investigate the influence of this attenu- 
ation on cosmological initial conditions in what follows. We 
analyze non-adaptive unigrid initial conditions in this section, 
in order to differentiate these effects from those due to adap- 
tive initial conditions, which will be addressed in Section 14.31 



3.4-1 Damping of small-scale perturbations 

In the Zel'dovich approximation, the displacement and veloc- 
ity fields are proportional to the gradient of the potential (cf. 

In this section, we investigate how the order of the 



11) 



eq. 

finite difference approximation for the Laplacian and the gra- 
dient affects perturbations close to the Nyquist wave number 
A^Ny = 7r//i, where h is the grid spacing which we set equal to 
unity in this section for convenience. 

We define v as the gradient of the potential arising from 
a source field /. The potential is solved using the mult i- grid 
scheme outlined above, and the gradient operator is applied 
subsequently. The exact solution v in k- space is given by v, 
where the tilde represents the Fourier transform. Expressed 
by the Fourier transforms of the operators, this is simply: 



(23) 



where the last equality holds for the exact solution. We define 
the attenuation as the ratio between the one-dimensional finite 
difference solution - using approximations of a given order for 
L and G and taking their Fourier transform as in Table [l] - 
and the exact solution v in k-space. 
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In Figure [S] we show this attenuation as a function of 
wavenumber k that is to be expected from the finite difference 
approximations to the operators up to sixth order. All finite 
difference gradients have zero amplitude at k^y so that fluctu- 
ations at this scale can not be represented in principle. A 2nd 
order approximation leads to significant attenuation of ^ 78 
per cent at k^y/2. However, for 6th order, attenuation is at 
the level of a few per cent at A:Ny/2. We thus expect a suppres- 
sion of the power spectrum at large k depending on the order 
of the finite difference operators employed. Note that k-^y/2 
corresponds to scales of two grid cells and that the attenuation 
enters squared into the power spectrum. 

3.5 Recovering small-scale power 

Both the use of a transfer function evaluated in real space and 
finite difference methods lead to a loss of power at scales close 
to the Nyquist wavenumber. We outline methods to solve this 
problem in what follows. Most of these methods will lead to 
spectral leakage and the associated spurious oscillations (cf. 
3.5.5) which poses no problem in the case of dark matter par- 



transforms of the finite difference operators are known, the 
fine grid solution can be simply replaced by the /c-space ex- 
act solution. Setting the right-hand side of Poisson's equation 
/(x) = on the boundary, we can recompute the self- gravity 
due to the fine grid and replace it with the one obtained with 
the exact /c-space method, i.e. using a grid zero-padded to 
twice its size, the correction is 



L(») 



/(k), 



(27) 



tide initial conditions but should be avoided for baryons. 



where n is the order of the finite difference approximation em- 
ployed and j is the direction along which the gradient is taken. 
The result is inverse transformed and added to the solution ob- 
tained with the finite difference method. The long-range part 
is still provided by the adaptive multi-grid solution, which is 
correct on scales larger than two grid cells. Thus, by defini- 
tion, small scale power is recovered, while the long-range part 
remains unaffected. Note that for the hybrid solver, the finite 
volume correction from Section |3.5.1| should not be applied 
since the /c-space Poisson solver is not a finite volume method. 
Instead, the method outlined in the next Section should be 
employed. 



3.5.1 Finite Volume Correction 



As demonstrated in Section |3.3| the multi-grid method is a 
finite volume approach. The solution is determined by com- 
puting flux balances across grid faces. This implies that the 
source field / is simply a cell average. Hence, implicitly, every 
discrete value Td (x) is a piecewise constant approximation to 
T7^(x) which is incorrect in the particle case. The piecewise 
constant averaging is given in terms of a kernel convolution 



^grid(xO = Ws{-X.) ★^cont.(x). 



(24) 



In order to restore the small-scale fluctuation amplitudes, 
it however suffices to perform a deconvolution with the cell 
assignment function - equivalent to nearest grid point (NGP) 
assignment, see e.g. iHockney & Eastwood (1981), 



Ws 



n ^ 



Xi + 



l-H{x^ 



(25) 



where H is the Heaviside step-function and h is the grid spac- 
ing. Since H has an algebraic form for its Fourier transform 



Ws{k)= n 



ki 



(26) 



we can perform this deconvolution in the Fourier domain and 



thus recover some of the sub-grid power - see alsojJing ( |2005 ) 
who use a similar procedure for power spectrum estimation. 
Note that - by virtue of restoring sub-grid power - the de- 
convolution leads to spectral leakage. The associated ringing 
will however be completely filtered out by subsequent finite 
difference operations. 



3.5.3 Averaging correction 

Even with the hybrid approach, a simple grid assignment of 
the real-space transfer function, determined as in Section [2. 2| 
will lead to damping at small scales since, unlike in the /c-space 
transfer function case, no sub-grid information is present. We 
can however restore this sub-grid power also for the hybrid 
case which itself corrects only the attenuation due to the finite 
volume approximation. 

We can compute each value T7^(x) from an average over 
sub-grid scales of the highest level i (imagined on the next 
higher refinement level 1 of the mesh) . The value on level £ is 
then an average over 8 cells at twice the resolution, equivalent 
to a convolution with the sub-grid kernel 



^sg(x)= n 



5i) \ Xi + - ] + Sb 



(28) 



where is the Dirac ^-function and h is the grid spacing. 
The discretized transfer function Toirijk) is thus given by av- 
eraging over the true transfer function T-jz at sub-grid scales 
Tnivijk) = Ttz Ksg. The kernel Ksg can be explicitly calcu- 
lated in Fourier space to be 



i^.g(k)=n< 



h ki 



(29) 



Power is thus reduced at large k compared to higher resolution. 
Deconvolving with Ksg will restore this power. The average 
over sub-grid cells is taken while computing the real space 
transfer function kernel on the three dimensional grid. 



3.5.2 A hybrid Poisson solver 

The attenuation of subgrid power due to the finite difference 
operators themselves (which however comes at the benefit of 
a non-oscillatory velocity field) may be considered undesir- 
able as part of the velocity information is effectively destroyed 
by the finite difference approximation. Since the lack of small 
scale power is only relevant on the finest grid, a simple solution 
to circumvent this problem can be devised. Since the Fourier 



3.5.4 Effect of the corrections on the initial power spectrum 

In Figure [6] we show the influence of the finite difference order 
on the initial power spectrum compared to traditional initial 
conditions using the k-space transfer function and a Fourier 
based Poisson solver. All initial conditions are for a 512"^ par- 
ticle discretisation of a 100/i~^Mpc box. The power spectrum 
was computed using the FFT on a 1024^ cell density grid ob- 
tained via cloud-in-cell (CIC) interpolation from the particle 
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Figure 6. Influence of the real space transfer function, the flnite 
difference approximations of the Laplacian and gradient operators 
and the sub-grid corrections on the initial power spectrum. Shown 
is the FFT estimated power spectrum (using CIC particle inter- 
polation on a 1024^ grid for the 512^ particles in a 100 h~^Mpc 
box) for the classical exact k-space initial conditions (solid black), 
uncorrected real-space, 6th order (dash-dot-dot) and the corrected 
version thereof (dash-dot), as well as the uncorrected hybrid (dot- 
ted) and the corrected hybrid (dashed). The solid gray line indicates 
the theoretical input spectrum. 



positions. No deconvolution with the grid assignment operator 
was performed so that the CIC assignment leads to attenua- 
tion close to the Nyquist frequency of the mesh (outside of the 
plot) used to compute power spectrum in all cases ( Jing|2005 ). 
The uncorrected spectra using the real-space transfer function 
kernel show lack of power for wave numbers above ^ 0.1 /cNy 
even when using the hybrid Poisson solver which is equivalent 
to a pure FFT Poisson solver in the case of a periodic grid. 
This lack of power on such large scales is undesirable as it will 
suppress the growth of well resolved haloes. 

Performing the finite volume correction in the case of the 
6th order finite difference Poisson solver recovers all power 
below ^ 0.5 /cNy but (naturally) retains the lack of power on 
even smaller scales. In the case of the hybrid Poisson solver, 
the power is recovered on all scales down to k^y. We observe 
a slight excess of power on the smallest scales. 

Note that the motivation to apply these corrections is to 
match the properties of /c-space sampled initial conditions as 
closely as possible. However, they should not be understood as 
mandatory and are treated as options in our code. Applying 
the corrections for the respective Poisson solvers, the initial 
power spectra can be made to agree with the traditional exact 
A:-space sampled initial conditions with differences only on the 
smallest scales. We investigate the influence of these differ- 
ences on A/'-body simulations in Section [4] For the remainder 
of the paper, the respective corrections are applied for each 
method without explicit reference. 
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Figure 7. Spatial distribution of the relative difference between 
the real space computed y-velocity fleld due to a single peak and 
the velocity fleld computed using /c-space sampling. The left panel 
shows the difference for the 6th order flnite difference fleld and the 
right panel for velocities determined with the hybrid solver. Use of 
the FFT poisson solver (in both the /c-space sampled and the hybrid 
case) leads to pronounced oscillations along the y-axis through the 
peak. The use of the /c-space transfer function leads to a loss of 
isotropy. Figure [s] shows the error along this axis in more detail. 
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Figure 8. Slice through the central part of the velocity fleld in- 
duced by a single peak (cf. Figure |7|. The top panel shows the 
y-velocity component using the three methods, the top panel shows 
the relative difference of the real-space methods with respect to the 
fc-space sampling. The left half of the top panel has positive sign, 
the right half is negative. 



3.5.5 Spectral leakage and aliasing in FFT-Poisson solvers 
and Fourier space transfer functions 

We investigate now the influence of using a real-space transfer 
function as well as a finite difference or hybrid Poisson solver 
compared to the standard method of using the /c-space transfer 
function and computing velocities. Velocities computed with 
equation (23) are not periodic in /c-space. 



Furthermore, the 
real space transform will be periodic and suffer from aliasing 
due to its Fourier space discretisation. 

In order to assess the associated effects, we compute the 
velocity perturbation associated with a single peak at the cen- 
ter of a 100/i~^Mpc box at 512"^ resolution. The peak repre- 
sents a Dirac ^-function convolved with the transfer function. 
In Figure |7| we show the difference between the velocity field 
obtained with the real space transfer function using both 6th 
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order multigrid and the hybrid Poisson solver, with respect to 
the velocity field obtained in the traditional way. We show two- 
dimensional slices in {x,y) of the relative difference through 
the peak for the y-component of the velocity. Note that the 
velocity field is zero on the horizontal axis through the centre. 

We first note that the relative difference falls below 10 ~^ 
only in the central region. The 6th order multi-grid result has a 
slightly smaller central region with such low difference to the k- 
space sampled results. Outside the central region, the relative 
difference rises to a maximum of ~ 2 per cent (outside the 
plane where the velocity vanishes), which is independent of the 
Poisson solver employed and thus a feature of the real-space 
transfer function. Indeed, it is the manifestation of the lack 
of spherical symmetry and aliasing (due to the finite number 
of samples) of the A;-space transfer function when transformed 
back to real space on a finite mesh. 

We also observe spectral leakage showing up as a vertical 
feature passing through the centre. Spectral leakage is due to 
the non-periodicity in Fourier-space and thus the truncation 
of a non-periodic function at a finite wavenumber. The trunca- 
tion leads to an oscillation at the Nyquist wavenumber in real 
space. In Figure [S] we show the central part of the velocity 
fields along the vertical axis of Figure [7| revealing the oscil- 
latory nature. Spectral leakage is present whenever a /c-space 
exact Poisson solver is employed. Hence, the feature is not 
present in the multi-grid solution, which is periodic in Fourier 
space. Naturally, a similar feature is also present in the /c-space 
transfer function itself, when transformed to real space on a 
three dimensional mesh (see also Bertschinger||2001 ). 

Note that leakage in the velocity perturbations associ- 
ated with a peak is likely to cause secondary correlations in 
the velocity field around rare high peaks but is unavoidable 
if no window function is used with the FFT Poisson solvers 
to truncate the spectrum smoothly at the Nyquist frequency. 
The use of such a window function (as e.g. the Hanning filter 
employed in Grafic-2, Bertschinger 2001) however reduces 
power at high frequencies which might be undesirable since it 
suppresses the small wavelength perturbations and thus the 
growth of haloes associated with them. We note that the fi- 
nite difference approach has similar filtering properties with 
a relatively sharp cut-off in /c-space. Furthermore the finite 
difference operators are periodic in Fourier space. We do not 
consider the use of other window functions but perform our 
analysis of initial conditions generated with finite difference 
multi-grid and the hybrid Poisson solver in parallel through- 
out the remainder of this paper. 




Figure 9. Redshift evolution of the power spectrum for initial 
conditions computed using a 6th order multi-grid Poisson solver 
(solid black line in top panel) and using a hybrid Poisson solver 
(solid black line in bottom panel) compared to the traditional k- 
space sampling approach (dashed red lines) and the power spec- 
trum evolution predicted by linear perturbation theory for the input 
power spectrum (gray lines). Results are shown for a simulation of 
a 1 h~^Gpc box with 512"^ particles. Power spectra have been com- 
puted after CIC mass assignment on a 1024^ grid. 



4.1 Statistical properties of the cosmological 
simulations 



4 ERROR ANALYSIS 

In this section, we first investigate the differences between the 
new methods described in Section [2] and |3] and the traditional 
approach of operating exclusively in Fourier space on various 
observables of cosmological A/'-body simulations evolved over 
cosmic time as well as of a single galaxy cluster in order to en- 
sure their proper performance in unigrid simulations. We then 
investigate the errors arising in nested "zoom-in" initial condi- 
tions generated with our method in both the initial conditions 
and in a re-simulation of a galaxy cluster. 



In order to assess differences between initial conditions gen- 
erated with the various methods, we ran Gadget-2 N-body 
simulation of a l/i~^Gpc and a 100/i~^Mpc box with 512^ 
particles from redshift z = 45 to z = 0. First order initial con- 
ditions (ILPT) were generated using the same random seeds. 
The power spectrum from which the initial conditions were 
generated in all cases has density parameters Qrn = 0.276 
and Qa = 0.724, a Hubble constant Ho = 70.3 km s~^ Mpc~\ 
a normalization ag — 0.811 as well as a spectral index of 
ris — 0.961 and was computed using the fitting formula of 



Eisenstein & Hu (1998). We use Plummer equivalent soft- 



ening lengths e of 0.09/i~^Mpc for the l/i~^Gpc runs and 
0.009 /i~^Mpc for the 100/i~^Mpc runs, respectively. These 
box sizes were chosen in order to probe both the highly non- 
linear regime, with the 100/i~^Mpc box, and the mildly non- 
linear regime with the 1 /i~^Gpc box. 
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Figure 10. The effect of the Poisson solver used to compute the 
initial conditions on the friends-of-friends halo abundance in N-body 
simulations at redshift z = 0. Shown is the number of haloes as a 
function of their mass and the number of particles they consist of 
in the top panels. The bottom panels show the relative difference 
with respect to the "exact k-space" method. The left panels are 
for a l/i~-'^Gpc simulation box, the right panels for a 100 h~^Mpc 
box. We show the results for a multi-grid Poisson solver using 2nd 
(dotted green), 4th (dashed blue), and 6th (dash-dotted red) order 
finite difference approximations, as well as for the hybrid Poisson 
solver (dash-dot-dot black lines) and the fc-space sampled initial 
conditions (solid gray). Errorbars in the lower panel indicate the 
Poisson errors on the halo counts. 



4.1.1 Evolution of the power spectrum 

In Figure [5] we show the redshift evolution of the density 
power spectra for the l/i~^Gpc box. We compare the linear 
theory evolution with the evolving power spectrum from ini- 
tial conditions generated with the traditional /c-space sampling 
method, as well as generated with a real space transfer func- 
tion convolution and both a 6th order finite difference and the 
hybrid Poisson solver. 

As previously discussed in Section [3. 5. 4| the finite differ- 
ence initial conditions show a significant suppression of power 
on the smallest scales. We see however that the power is im- 
mediately restored once the relevant scales enter non-linear 
growth and power is transferred from larger to smaller scales. 
The difference with respect to the /c-space initial conditions 
decreases over time and is no longer visible at ^ ~ 1. No 
wavenumbers below ~ 0.5A:Ny of the 512^ initial particle grid 
are affected at any time. Initial conditions generated with the 
hybrid scheme agree (almost) perfectly with the /c-space ver- 
sion. Again, the small initial difference also disappears com- 
pletely once the relevant scales become non-linear. 

Regarding large scales, we see no difference between the 
three methods, so that we focus our attention on a closer in- 
spection of small scale behaviour in what follows. 



4.1.2 The z = halo abundance 

We compare halo abundance at a given mass between the var- 
ious initial conditions. Haloes are identified using the Friends- 
grouping 



of-Friends (FoF) algorithm ( [Press k Davis] |1982| ), 
particles which are separated by less than 0.2 times the mean 
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Figure 11. Same as Figure^Ojbut showing the effect of the Poisson 
solvers used to compute the initial conditions on the small-scale two- 
point correlation function at redshift z = 0. Errorbars in the lower 
panel indicate the Poisson errors. 



inter-particle distance. Figure shows the number of haloes 
A^haioes iH the simulatious as a function of both halo particle 
number A^part and halo mass. For both boxes and as expected, 
we see the strongest difference between the /c-space sampled 
initial conditions and those with 2nd order accuracy of the 
finite difference operators. In fact, halo abundance is lower up 
to 1000 halo particles. The histograms agree above ^100 and 
above ^ 40 halo particles for 4th and 6th order, respectively, 
with that of the /c-space method. The hybrid initial conditions 
agree over all mass ranges within ^ 2 per cent. Since halo 
abundance deviates from the Schechter function form also for 
the k-space initial conditions for haloes below 100 particles, it 
remains a question of future investigation whether the over- 
abundance of the smallest mass haloes is a spurious result of 
the /c-space initial conditions. 



4.1.3 The z = small-scale correlation function 

In Figure [11] we show the two-point correlation functions at 
z = 0. The two-point correlation function ^(r) is defined as 
the excess probability (with respect to random) to find two 
particles in volume elements dVi and dV2 separated by a co- 
moving distance r. 



dPi2(r) = [l + eWldFidys, 



(30) 



where n indicates the mean particle number density. We com- 
puted the two-point correlation function using tree based 
sparse sampling thus computing more close-pairs than distant 
pairs for a subset of 1 per cent of the particles. Thus, statistical 
errors increase on larger scales. 

For the l/i~^Gpc box, our results indicate that ^ is re- 
duced by ^ 17 per cent at r = e when using a 2nd order finite 
difference scheme, while for 6th order the reduction is only 
around ^ 5 per cent. Both 4th and 6th order results deviate 
by less than ^ 2 per cent at scales r > 6e. For the highly 
non-linear regime probed by the 100/i~^Mpc runs, we do not 
find any significant differences. Gravitational interaction has 
wiped out any initial difference. We conclude that differences 
are negligible in the highly non-linear regime and as soon as 
mildly non- linear scales are well resolved. 
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Figure 12. Galaxy cluster of mass 2 x 10-*^^ h~-^MQ at z = for the 
different methods to set up initial conditions. The top panel shows 
the cluster formed from traditional fc-space sampled initial condi- 
tions, while the bottom panels correspond to initial conditions using 
a convolution with the real space transfer function in conjunction 
with a 6th order finite difference Poisson solver (bottom left) and 
the hybrid Poisson solver (bottom right). Each panel is the projec- 
tion of a cube of 2 x 2 x 2 h~^Mpc centered on the densest particle 
in the cluster. 



4.2 Analysis of a. z — galaxy cluster halo 

In order to assess the influence of the real-space approach on 
the radial density profile of haloes, we investigate in more de- 
tail the properties of one cluster halo forming from initial con- 
ditions generated in the various ways in the 100/i~^Mpc box 
with 512"^ particles. 

In Figure [12] images of the cluster at z = are shown 
for the three methods. Contrast in the images has been set in 
such a way as to make substructure most visible, white areas 
are not actually devoid of particles. Visually, the only differ- 
ence between the three panels lies in the position of smaller 
substructure within the halo. All larger sub-haloes agree in 
their position. The overall shape of the halo is identical. 

In order to asses differences more quantitatively, we show 
in Figure [13] the radial density profile centered on the dens- 
est particle (determined by averaging over 32 neighbours and 
weighted by a spline kernel) in each case. We find no system- 
atic difference between either of the approaches. Differences 
are around 5 per cent and show no obvious bias. They are 
likely attributable to slight differences in the position of the 
substructure in the halo. These results further support the 
conclusion made in the previous section that differences be- 
tween the real and the k-space approach are completely neg- 
ligible in the highly non-linear regime. 

The gross properties of the cluster halo are analyzed using 
the Amiga halo finder (AHF) ( jKnollmann k KnebepOOOt . 
Table [3] we give the virial mass Mvir and radius i^vir, max- 
imum circular velocity Vmax, the radius Rmax at which the 
circular velocity profile is maximal, the dimensionless spin pa- 
rameter A ( Bullock et al.|20"oT ), the three-dimensional velocity 
dispersion ay as well as the minor-to-major axis ratio c/a, de- 
termined using the eigenvalues of the moment of inertia tensor 
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Figure 13. Radial density profiles for a 2 x IQ-*^^ h~^MQ halo at z = 
0. The top panel shows the profiles for initial conditions obtained 
with the 6th order multi-grid (dotted blue) and the hybrid Poisson 
solver (dashed green), as well as for the /c-space sampled initial 
conditions (solid black line). The bottom panel shows the relative 
difference between the profiles from real space and k-space initial 
conditions for the two methods. 



T(/c), fc-space T(r), finite, diff. T(r), hybrid 



Mvir//l-^M0 


2.140 X IQi^ 


2.141 X 10^4 


2.130 X 10 


^vir / h~^'kpC 


1231.2 


1231.4 


1229.4 


Vmax /kmS~^ 


941.5 


929.4 


940.6 


i^max / h~^kpC 


680.6 


648.0 


669.5 


X 


0.02510 


0.02497 


0.02503 


CTv / kms"-*^ 


1000.1 


986.0 


1001.5 


c/a 


0.7119 


0.7119 


0.7144 



Table 3. Properties of a cluster halo at z = for different Pois- 
son solvers used to compute the initial conditions and real//c-space 
sampled density fields. 



(see e.g. Hahn et al.||2007 ). With the exception of i?max, dif- 
ferences are around 1 per cent for all quantities investigated. 
It appears that i^max is very sensitive to small changes in the 
density profiles leading to a difference of a few per cent in the 
finite difference case. 

Finally, using AHF, we investigate the sub- halo popula- 
tion of the cluster. Since the finite difference approach sup- 
presses slightly the formation of the smallest haloes (cf. 4.1.2), 



it is plausible that sub-haloes are affected as well, possibly even 
stronger. In Figure [14] we show the cumulative distribution of 
maximum circular velocity t'max for all sub-haloes with at least 
30 particles contained within the virial radius of the main clus- 
ter halo. We observe that the most massive substructure has 
a slightly lower Vma^ for the /c-space initial conditions than for 
the others. The following sub-haloes, at next lower Vma^ are 
identical. For the lower mass sub-haloes, the cumulative dis- 
tributions agree, the simulation starting from the mult i- grid 
initial conditions having slightly larger scatter than the one 
starting from hybrid initial conditions. In particular, we see 
no systematic deficiency of lowest mass haloes in the multi- 
grid case as might have been expected. 

We conclude that differences in halo properties between 
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Figure 14. Cumulative distribution of maximum circular velocities 
Vmax of sub-haloes of the cluster for the various methods. 

the real-space and /c-space conditions are at the per cent level 
for gross properties and at the few per cent level for profiles 
and sub-halo mass functions. We do not observe any system- 
atic bias. 



4.3 Errors in nested initial conditions: Initial 
redshift 

In this section, we analyze the errors in nested initial condi- 
tions computed with our method compared to a full-resolution 
solution with no local refinement. In order to compare our re- 



sults with Bertschinger (2001), we decided to implement an 



identical set-up, where we replace white noise in cells outside 
the refinement region in the full-resolution set-up with the re- 
spective average value it has on the coarse grid in the "zoom- 
in" set-up. This ensures that identical white noise information 
is used in all cases, i.e. all the information that is in princi- 
ple representable in the refined hierarchy. Note that naturally 
errors must arise at coarse-fine boundaries since the velocity 
and displacement fields must transition smoothly from fine to 
coarse resolution. 



4.3.1 One refinement level 

First, initial conditions are generated on a uniform 512^ grid 
for a 100/i~^Mpc box with white noise degraded to include 
identical information as in the "zoom-in" set-up to obtain 
the full-resolution answer. Next, another set of initial con- 
ditions is generated for a region of 0.2 times the box length 
at an effective resolution of 512^ while the remaining volume 
is treated at a resolution of 256^. The high-resolution region 
grid cells. In Figure 



15 



we show the error for 



comprises 102 

the density field and the three components of the velocity for 
a slice in the x — y plane through the centre of the refinement 
region. Colors correspond to the logarithm of the difference 
between the multi-scale solution and the full-grid solution in 
units of the standard deviation erg of the respective quantity 
Q ^ {S,Vx,Vy,Vz} on the full grid. In order to suppress the 
contribution from density convolution errors at the bound- 
ary, we padded the grids with additional 8 grid cells at each 
boundary during the convolutions. This padding is however 
cut away before the potential is computed for the multi-scale 



hierarchy to ensure that the final velocity field is continuous 
and different iable across coarse-fine boundaries. All results are 
shown for ILPT initial conditions, the respective errors in the 
2LPT case are qualitatively identical. We used the 6th order 
multi-grid for the results shown in this and the next section. 

For the density, almost everywhere the difference is sig- 
nificantly below lO^^cr^. We see that some errors at a few 
^ 10~^as resulting from the interpolation of long-range com- 
ponents onto the finer grid. We find an RMS error of ^ 
7 X lO-^cr^. 

For the velocity components, again, errors are signifi- 
cantly below 10~^(7^;. in the interior of the refinement region, 
apart from a few cells surrounding the boundary. Errors at the 
boundary are at most ^ 10~^" ^cr^. and arise due to the smooth 
coarse-fine transition in the refined set-up which is not present 
in the unigrid case. The RMS error is ~ 1.5 x 10" V^. for all 
three velocity components. Since the interior errors are be- 
low 10"^cr, the RMS error over the entire refinement region is 
clearly dominated by only the errors at the boundary. 

Errors for the hybrid Poisson solver are harder to com- 
pare since the long-range component from coarse grids is taken 
from the finite difference multi-grid solution so that it has no 



spectral leakage (see Section 3.5.5), while the same long-range 



component in the unigrid situation will show the spectral leak- 
age in the form of oscillations at the Nyquist frequency around 
the smooth solution. This lack of spectral leakage in the long- 
range components completely dominates the error in the inte- 
rior in a direct comparison of the unigrid and the multi-scale 
solution at a level of ^ lO^^cr^^, no other errors are however 
introduced. Since leakage is an artifact of the finite mesh used 
for the FFTs, these differences are spurious and of no signifi- 
cance to our multi-scale approach to initial conditions. 

In Figure |16| we show, for comparison, the error in the 
density and Vx velocity component obtained with Grafic2 for 
the identical refinement set-up and identical random numbers 
as in Figure [15] We have used version 2.101 of the code and not 
removed or modified the standard antialiasing filter. Errors in 
density are about one, errors in velocity (and thus displace- 
ments) are about two orders of magnitude larger. Furthermore, 
unlike with our approach, errors are oscillatory and spread 
over the entire refinement domain. In addition, there is a sig- 
nificant vertical gradient at the level of several percent. We 
can only speculate about the sources of these errors. For the 
density field, our approach avoids the problem of aliasing and 
spectral leakage completely by operating purely in real-space. 
Furthermore, it produces by definition identical results for the 
^feif component in the refined and the unigrid case. In addi- 
tion, we have ensured that mass is conserved between grids 
thus avoiding long-range errors in the velocity components. 
Finally, the use of the adaptive multi-grid Poisson solver to 
compute velocities guarantees a solution of Poisson's equation 
that is independent of the presence of a refined region if the 
mass field in the unigrid case is averaged with the restriction 
operator and then re-injected with the respective injection op- 
erator (averaging and straight injection in our case) in order 
to refiect the resolution in the refinement hierarchy. 

As mentioned before, Grafic2 also uses an antialiasing 
filter which severely suppresses power on small scales, so that 
the density and velocity fields generated are very smooth. We 
conclude that our new approach is an improvement by about 
two orders of magnitude in error reduction and the additional 
advantage that errors are completely confined to the bound- 
aries. 
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Figure 15. Error for the over-density (left-most panel) and the three components of the velocity (vx, Vy, Vz - second to fourth panel, 
respectively) arising from the nested multi-scale convolution approach for one additional refinement region. Shown is the difference between 
the respective quantity computed using nested convolution on a 0.2 x 0.2 x 0.2 = 102^ sub- volume at an effective resolution of 512^ embedded 
in a coarse grid of 256^ cells and the quantity in that region computed using the full grid of 512^ cells. We show the base 10 logarithm of the 
difference in units of the standard deviation of the respective field inside the refinement region on the 512^ grid. The slices show the x — y 
plane through the centre of the refinement region. 




Figure 16. Error for the density and the Vx velocity component us- 
ing Grafic2 for the same refinement set-up as in Figure [Ts] Again, 
the base 10 logarithm of the difference in units of the standard de- 
viation of the respective field inside the refinement region is shown. 
Note that the colour scale for the Vx error corresponds to a range 
one order in magnitude larger than in Figure [Ts] The errors in the 
other velocity components are comparable to that for Vx- 

In summary, with our method, errors in the interior are 
below l^~^a with larger errors occurring at the boundary for 
the velocities due to the coarse- fine transition in the refinement 
hierarchy. It is particularly interesting to note that the velocity 
errors have no visible small scale components in the interior. 

4.3.2 Two refinement levels 

In this Section, we repeat the error analysis from above for a 
set-up of two nested refinement levels. The refinement region 
still has 102"^ cells at an effective resolution of 512"^ but is now 
embedded in an intermediate grid at 256^ effective resolution 
which is in turn embedded in the coarse grid with 128*^ cells 
extending over the entire volume. The intermediate grid has 



an extent of 16 cells beyond the edges of the finest grid. As be- 
fore, we pad each refinement level with additional 8 boundary 
cells during the convolution step in order to avoid coarse-fine 
transition errors during the convolution. These boundary cells 
are again cut off before the potential is computed. 

In Figure |17[ we show the differences in the high- 
resolution region between the full-grid solution and the nested 
multi-scale solution. Density and velocity errors in the high- 
resolution region are roughly identical to the one level case 
above and below lO"'^ ctq in the interior. We find RMS errors 
of ~ 5 X 10~^cr(5 for the overdensity field and 1.8 x 10""^ cr^^ 
for the three velocity components. 

To summarize, we see no significant impact on the veloc- 
ity field due to the introduction of additional refinement levels. 
All errors are completely localized at the coarse-fine bound- 
aries. Regarding memory consumption, using double preci- 
sion variables (8 bytes), we observed a peak memory usage of 
~ 2.7 GBytes in unigrid mode (a 512^ double precision array 
requires 1 GByte of memory), ^ 670 MBytes in the case of the 
single refinement level set-up from above and ^ 270 MBytes 
in the two-level case from above, demonstrating the low mem- 
ory footprint of our "zoom-in" approach to generate initial 
conditions. 



4.4 Errors in nested initial conditions: 

Re-simulation of a galaxy cluster halo 

In order to assess the influence of multi-scale initial conditions 
on the formation of a re-simulated object, we generate multi- 
scale initial conditions for one of the most massive clusters 
with mass 2 x 10^^ /i~^Mpc in the 100/i~^Mpc box discussed 
before. We deliberately set the rectangular refinement region 
to include only the Lagrangian patch of the cluster halo. This 
will maximize the influence of both the coarse sampling of the 
large-scale tidal field and boundary effects on the formation of 
the cluster halo and allow us to estimate these effects. The re- 
finement volume is a rectangular region of 23 x 21 x 21 /i~^Mpc^ 
and was determined by following the FoF particles consti- 
tuting the halo at z = back to the initial conditions and 
determining their bounding box rounded up/down to integer 

Mpc. Note that a sphere at mean density containing the 
mass of the cluster would have a diameter of ^ 17.5/i~^Mpc. 

We perform Gadget-2 simulations with one level of re- 
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Figure 17. Same as Figure [is] but with two refinement levels. The 102^ grid at 512^ effective resolution is embedded in an intermediate grid 
at 256"^ effective resolution which is again embedded in a coarse grid with 128"^ cells. The intermediate grid pads the fine grid with 16 cells 
on each side. The errors are shown only for the finest grid. 
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Figure 18. The resimulated cluster at 2; = with 512^ effective resolution. Panels in the upper row were obtained by degrading the initial 
conditions outside the refinement region, while initial conditions for the bottom row panels were generated using the multi-scale scheme. Left 
panels use the 6*^ order finite difference scheme, right panels the hybrid Poisson solver. In each of the two panels, the left column corresponds 
to a simulation set-up with 256^ effective resolution outside the high resolution region, while the right column shows result for a two-level 
zoom with a bounding region at 256^ resolution around the high resolution region and 128^ resolution in the remainder of the box. Each 
image represents a projected cube of 2 x 2 x 2/i~^Mpc^, only high-resolution particles are shown. 



finement where the effective resolution in the high resolution 
region is 512*^ and 256"^ in the remainder of the box, as well 
as with two levels of refinement, where the high resolution re- 
gion is surrounded by a layer of 16 particles thickness around 
each face at 256"^ effective resolution that then drops to 128^ 
in the remainder of the box. In addition, we perform each sim- 
ulation two times for both the multigrid finite-difference and 
the hybrid approach: once with the initial density field deter- 
mined using the multi-scale convolution technique described 
in Section [23) and once with the density field determined at 
full resolution and then degraded (by averaging) to the same 
resolution as in the multi-scale setup followed by solving Pois- 
son's equation on the nested grid rather than the full grid. 
The latter produces completely negligible errors in the ve- 
locities and displacements inside the high resolution region. 
Differences from the full-grid initial velocities occur only in 
2-3 cells at the boundary where the fields transition smoothly 
from fine to coarse resolution. All simulations are run with 
Gadget-2 and use a force softening length of 0.009 /i~^Mpc 



for the high-resolution particles and 0.09 h ^Mpc for the other 
particles. 

In Figure [18] we show the clusters ai z = for all combi- 
nations of refinement set-up, Poisson solver and initial density 
generation method. We observe that the main halo as well as 
the most massive subhalos are consistent in position and size. 
The positions of some of the smaller halos are shifted and some 
smaller haloes seem to have merged in one set-up while they 
have not in some of the others. In general, visual differences 
are minimal. It is furthermore surprising that no systematic 
difference between the 1-level and the 2-level set-up can be 
seen. This is most likely a result of the padding region in the 
2-level simulations at 256^ resolution which is identical to the 
1-level simulations - the tidal influence of structures outside 
the padding region thus appears negligible and highlights the 
importance of adding padding. 

Comparing Figure ^] with Figure |12[ which shows the 
cluster in the unigrid simulations, we observe no obvious sys- 
tematic differences apart from smaller sub-halo positions be- 
ing slightly shifted. The overall shape of the cluster halo agrees 
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Figure 19. Overdensity profiles for the cluster at z = re- 
simulated with one and two refinement levels (top). The lower panel 
shows the relative difference with respect to the profile of the cluster 
in the unigrid simulation. Left panels show the result obtained from 
finite difference multi-grid initial conditions, the right panels corre- 
spond to the hybrid approach. The vertical light gray lines indicate 
the scale of three times the softening length. 
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Figure 20. Cumulative histogram of sub-halo maximum circular 
velocities Vma^ of the re-simulated cluster at 2; = 0. The left panel 
shows the result for finite difference multi-grid approach, the right 
panel the corresponding results for the hybrid approach. 



very well between the re-simulations and the full box simula- 
tions. 



We ran the Amiga halo finder (AHF) ( Knollmann 
Knebe|[2009 ) on the resimulated clusters and quote the same 
key values as in Section [4j] to quantify the gross properties of 
the cluster halo in Table [4] As in the unigrid case, apart from 
i^max, we find differences around 1 per cent or below for all 
quantities investigated. 

In Figure [19] we show the radial over-density profiles for 
the various cases. The lower panel show the relative difference 
with respect to the unigrid simulations. We observe no bias in 
either case. Scatter around the unigrid profile is larger for the 
finite difference case than for the hybrid initial conditions. For 
both methods, the scatter is slightly larger for the 2-level than 
for the 1-level set-up. Note that all "zoom- in" density profiles 
fall below the unigrid profiles in the last bin. This is due to the 
too small size of our refinement region as it is also present in 
the simulations generated by degrading the full density field. 
Despite the fact that it had been chosen too small, we find 
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Figure 21. Radial density profiles of the re-simulated cluster haloes 
at increasing resolution with a base grid of 128"^ particles. The pro- 
file for 1 refinement level, 256^ effective resolution is shown as a 
dotted black line; 2 levels, 512^ effective resolution as a dashed 
red line; and 3 levels, 1024^ effective resolution, as a solid blue line. 
The vertical light gray lines in the corresponding line styles indicate 
three times the softening length for each simulation. 



however not a single low resolution particle inside the virial 
radius, the first appearing at ^ 1.3i?vir- 

In Figure [20l we show the abundances of substructure as 
a function of maximum circular velocity t'max within the virial 
radius of the cluster for all of the re-simulations. Differences 
between the 1 and 2-level results are at the level of several per 
cent for the finite difference case. The hybrid initial conditions 
show better agreement between the four runs. It is hard and 
beyond the scope of this paper to investigate to what degree 
these differences stem from simple changes of the sub-halo 
positions causing errors in determining their circular velocities 
(or masses) in the sub-halo finder. 

We observe that the errors due to our multi-scale method, 
particularly with the hybrid Poisson solver, agree very well 
with the degraded initial conditions in both the one and two 
refinement level set-up. In particular, the scatter between the 
results from multi-scale density fields and degraded density 
fields is never larger than the difference between degraded ini- 
tial conditions at one or two refinement levels itself. This leads 
us to conclude that the differences, apart from the introduction 
of an additional stochastic component, are dominated by the 
late-time evolution and not the initial conditions. In particu- 
lar, we find no evidence that our method introduces systematic 
difference or bias in any of the investigated quantities. 

Finally, we investigate the convergence of the density pro- 
file at even higher resolution for the hybrid Poisson solver case. 
In Figure |21| we show the radial density profiles for a series 
of re-simulations that all have a base resolution of 128^ parti- 
cles and one to three additional refinement levels. The profiles 
trace each other almost perfectly. 
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Mvir //^~^M0 


Rvir / h ^kpc Vmax/kmS ^ Rma^ / h ^kpC 


A 


(Tv / km s ^ 


c/a 


6th order finite difference 


1 level, 


degraded 


2.136 X IQi^ 


1230.4 


939.5 


646.5 


0.02543 


998.5 


0.7168 


1 level, 


multi-scale 


2.137 X 10^4 


1230.6 


941.4 


667.7 


0.02545 


999.3 


0.7168 


2 level. 


degraded 


2.138 X 10^4 


1230.9 


938.5 


652.3 


0.02552 


997.4 


0.7141 


2 level. 


multi-scale 


2.137 X 10^4 


1230.6 


939.8 


612.8 


0.02550 


998.3 


0.7171 


hybrid Poisson solver 


1 level. 


degraded 


2.135 X 10^^ 


1230.3 


942.4 


667.9 


0.02525 


1000.7 


0.7142 


1 level. 


multi-scale 


2.134 X 10^4 


1230.1 


941.4 


649.8 


0.02530 


999.9 


0.7136 


2 level. 


degraded 


2.134 X 10^^ 


1230.0 


940.2 


665.1 


0.02532 


998.4 


0.7130 


2 level. 


multi-scale 


2.134 X 10^4 


1230.0 


941.8 


669.6 


0.02534 


999.5 


0.7145 



Table 4. Properties of the cluster halo when resimulated with one or two levels of coarse resolution particles outside the high-resolution 
region. The 'degraded' runs started from initial conditions that were computed at the full resolution of the high-resolution region, while the 
'multi-scale' runs use the adaptive zoom-in technique. 



5 BARYON INITIAL CONDITIONS 

In this section, we discuss the generation of initial conditions 
for a two- component - dark matter and baryon - fluid. Linear 
perturbation theory predicts distinct amplitudes for baryon 
and CDM density fluctuations and also for their respective 
velocity fluctuations. We demonstrate that our approach using 
real-space transfer functions can be easily generalized for such 
a two-component fluid. 



5.1 Density and velocity transfer functions 

At high redshifts, at which initial conditions for cosmological 
simulations are typically generated, baryon density fluctua- 
tions are not yet exactly tracing the dark matter fluctuations 



tions. We do not include the baryon temperature fluctuations 



(e.g. Yamamoto et al. 1998). Furthermore, baryon fluctua- 



tions are exponentially suppressed beyond the photon diffu- 
sion scale ( Silk| p"968 ) . These effects lead to a different shape 
of the baryon transfer function compared to the pure dark 
matter transfer function, particularly on small scales below 
- 10"^/i"^Mpc. 



[Yosh ida et al.| ( 2003 ) demonstrated that the growth of 
density perturbations in the two-component fluid can only be 
correctly reproduced if besides the different initial amplitudes 
of density perturbations also the difference in initial velocities 
between the two components are respected. It is thus impor- 
tant that initial conditions for the two-component fluid reflect 
these important differences between baryons and dark matter 
and are thus able to reproduce the correct growth of fluctua- 
tions in both components consistent with the predictions from 
linear perturbation theories for those scales where perturba- 
tions are still small also with cosmological A/'-body + hydrody- 
namics simulation codes. We defer a more detailed comparison 
with non-linear perturbation theory in two-component fluids 
(e.g. Somogyi Smitli||2010 ) to later work. 

For the two component fluid of baryons and dark mat- 
ter, we determine both baryon and dark matter velocities us- 
ing transfer functions, equivalent to the density perturbation 
transfer functions, that we construct in the following way. Us- 
ing a Boltzmann solver, we solve the linearized equations of 
relativistic perturbation theory to obtain the velocity pertur- 
bations of the CDM and the baryon fluid. We adopt the nota- 



tion of Ma & Bertschinger ( 1995 ) and kindly refer the reader 
to that publication for all details on the perturbation equa- 



as suggested by Naoz et al. (2010) as it is mainly relevant to 
perturbations with wavenumbers at and above the strongly 
Jeans damped regime. Using the results of linear perturbation 
theory, we define the two velocity transfer functions 



— QjuHq 



a 



(31) 
(32) 



where (j) is the conformal Newtonian potential for spatial met- 
ric perturbations and Ob is the baryon velocity divergence. 
These transfer functions are constructed in such a way that at 
linear order the relation 

i;,(fc) oc ^T4fc)r»/^ i = 1,2,3 (33) 

holds for the amplitude of velocity perturbations, i.e. they 
simply replace the density transfer functions in our algorithm 
when computing initial velocities as in Section [3T7] They 



furthermore allow the application of 2LPT (cf. Section 3.1.2), 
as they generate an "effective" source field for the velocity 
perturbations completely analogously to the relation between 
density and displacement. The respective real space transfer 
functions for all four fields are shown in Figure [22] 



5.2 Evolution of Baryon perturbations 

We now analyze the evolution of distinct baryon and CDM 
perturbations in two-component simulations using three 
commonly used cosmological simulation codes: Gadget-2 
(Springel 200 5|), Ramses (|TbyssTer|2002|) and Enzo ([Bryan 
,Norman„1997f Bryan et al.[[2001[ [O'Shea et al.[[2004t . In each 



case, the initial conditions are generated at 512 resolution for 
a box of comoving size 8/i~^Mpc with the same cosmological 
parameters as before and r^baryon = 0.045. The initial redshift 
is set to 2; = 100. For the grid codes (Enzo and Ramses), we 
do not allow for adaptive mesh refinement in order to avoid 
multi-scale density fields when computing the matter power 
spectra. In the SPH case (Gadet-2), we use a gravitational 
force softening of comoving 0.8/i~^kpc for both baryons and 
CDM. Furthermore, SPH particles are placed on a staggered 
grid with respect to the CDM particles. In all cases, we as- 
sume a polytropic equation of state with adiabatic exponent 
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Figure 23. Evolution of the density power spectra of CDM (solid) and baryon (dashed) perturbations for different cosmological simulation 
codes. The initial conditions have distinct amplitudes for baryon and CDM density and velocity perturbations. We show the power spectra 
at redshifts z = 100, 50, 30 and 20 (bottom to top), the corresponding results from linear perturbation theory are shown as gray lines in the 
background. Note that no adaptive mesh refinement was allowed for Ramses and Enzo. 
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Figure 22. The real space transfer function T(r) at z = 100 for 
baryon (red) and CDM (blue) density and baryon (green) and CDM 
(orange) velocity perturbations. The inset shows the first baryon 
oscillation peak in linear scale. Note the effect of the Jeans scale for 
r < 10~^ h~^Mpc on the baryons. 



7 = 5/3, an initial temperature of 140K//i, where /i is the 
mean molecular weight, and use neither cooling nor energy 
feedback. Note that Ramses, in our set-up, uses a piecewise 
linear method (PLM) for the hydrodynamics and mixed 2nd 
and 4th order gravity (2nd order in the Laplacian, 4th order 
in the gradients). Furthermore, we use the multidimensional 
monotonized central slope limiter ( [van Leer|1977| ). Enzo uses 
a piecewise parabolic method (PPM) for the hydro and 2nd 



order gravity. We stop the simulations at ^ = 20 and thus 
probe into the mildly non-linear regime. 

The power spectra of the two-component perturbations 
are shown in Figure [23] at the initial and three subsequent 
redshifts, evolved with the respective codes. These power spec- 
tra have been computed on a 512^ grid using FFT with CIC 
interpolation for all particles, and directly from the baryon 
overdensity field for the grid codes. We do not correct for the 
loss of power due to the interpolation scheme and also do not 
smooth the SPH particles with their respective kernels. Thus, 
the baryon perturbations in the SPH case do not reflect the 
density perturbations seen by the SPH scheme, resolution thus 
appears to be significantly higher than for the grid codes, so 
that our results should not be used as an argument in favour 
or disfavour of either method. 

We observe several differences between the codes. The 
SPH run with Gadget- 2 shows excellent agreement with lin- 
ear perturbation theory. The relative amplitude between CDM 
and baryon perturbations at any given mode matches very 
well even in the mildly non-linear regime at 2; = 20. It is not 
suprising that Gadget performs very well in our test since hy- 
drodynamical forces are negligible in the regime probed and 
we thus only observe the performance of the tree-PM gravity 
solver. 

In the case of Ramses and Enzo, CDM perturbations on 
the smallest scales grow slower than predicted by linear PT. 
This is however expected and a direct comparison with the 
Gadget results is unfair since we did not allow for refinement, 
and gravity forces are thus smoothed at the grid scale. For 
Enzo, the growth of CDM perturbations is more damped at 
small scales than in the Ramses run. This is most likely the 
result of the higher order of the gradient used by the Ramses 
Poisson solver. 

For the baryon perturbations, small scales are growing 
slightly slower with Ramses than with Enzo. Again, this is 
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to be expected since we use Enzo with PPM and Ramses 
with PLM. Also the type of Riemann solver and slope lim- 
iter used in Ramses has a slight effect, more diffusive limiters 
and solvers slightly decrease the power on the smallest scales. 
On large scales, all codes reproduce the correct growth of the 
perturbations, in very good agreement with linear PT. 

To summarize, we observe that all three codes correctly 
evolve a baryon-CDM two-component fluid which starts with 
distinct density and velocity perturbations for each compo- 
nent. Differences that we observe with respect to the results 
from linear perturbation theory are attributable to the vari- 
ous numerical methods used in the simulation codes and not 
to the initial conditions. 



5.3 Local Lagrangian approximation for the baryon 
density field 

The density field sourcing the displacements and velocities in 
Lagrangian perturbation theory - which is a Gaussian ran- 
dom field - is inconsistent with the density field of the dis- 
placed particles - which is non-Gaussian. Using the Gaussian 
perturbation field as the initial baryon density field in grid 
based codes is thus inconsistent with the CDM perturbations. 
In this section, we describe an approach to prescribe initial 
gas density for the mesh cells based on the local Lagrangian 
approximation (LLA) that leads to consistency between the 
initial gas and dark matter density fields (see also Betancort 



Rijo|1991||Padmanabhan Subramanian|1993| [Protogeros 
Scherrer|1997| . Furthermore, it provides a natural way to pre 
scribe the initial gas density on a grid also for 2LPT initial 
conditions. 

The continuity equation is trivially fulfilled for a La- 
grangian description of the fluid and simply reads dm = const. 
We can however express the continuity equation in terms of 
the evolving fluid coordinates x. It then becomes p(x, t) d^x = 
pd^q, where p is the unperturbed mean density and p{:xi,t) is 
the density of the fluid element. Since the continuity equation 
is simply a volume integral, we can employ the change of vari- 
ables theorem to express the left-hand-side also in terms of 
the initial position q. 



^x 
9q 



<fq = pd^q. 



(34) 



Then, using eq. (10), the formal evolution of the density at 



the initial position of the fluid elements can be simply written 
as 



p(q,*) 



det 



(35) 



where Sij is the unity matrix. This implies that p(q, t) assumes 
a non-Gaussian distribution in general. In fact, p(q, t) will be- 
come singular whenever an eigenvalue of dLi / dqj will become 
— 1, corresponding to shell crossing along the respective axis. 
Since we use a first or second order Lagrangian perturbation 
theory approximation for the displacement field L(q, t), eq. 
( 35 ) is not exact and amounts to a "local Lagrangian approxi- 
mation". This means that, in general, mass will not be strictly 
conserved. We correct this by enforcing mass conservation by 
an a posteriori renormalization p ^ p-\- P: where p is the mea- 
sured deviation from mean density. Typically, for sufficiently 
high initial redshift (as e.g. for the z — 100 initial conditions 
used), the relative error is around 0.6 per cent and thus neg- 




5/a 



Figure 24. Evolution of the probability distribution functions 
(PDFs) of the baryon overdensity fields in units of the RMS over- 
density. For the 8 h~^Mpc simulation box, the PDFs are given at 
three redshifts: z = 50 (top panel), 2; = 30 (middle panel) and 
z = 20 (bottom panel). Initial conditions were generated at two 
starting redshifts Zs = 100 (solid lines) and Zs = 200 (dotted lines) 
with first order Lagrangian perturbation theory for dark matter and 
an initial Gaussian baryon density field following Eulerian pertur- 
bation theory (EPT; red), with first order LPT for dark matter and 
baryons (using the LLA) (ILPT; blue) and with second order LPT 
for dark matter and baryons (2LPT; green). 



ligible (note that the error amounts to no more than ^ 2 per 
cent at much lower redshifts). 

When computing L(q, t) for baryons, $ is sourced only 
by the baryon density perturbations, and so also ^ contains 
only baryonic contributions. Note that application of the LLA 
does not change the baryon power spectrum. 

In Figure [24] we show the evolution of the probability 
distribution function (PDF) of the baryon overdensity field in 
an Enzo simulation using the same simulation parameters as 
above but for two different starting redshifts Zs = 100 and Zs = 
200. For both starting redshifts, we generate initial conditions 
in three different ways: 

EPT: First order Lagrangian perturbation theory is used to 
initialize dark matter displacements and velocities as well as 
baryon velocities. The initial baryon density field is taken from 
Eulerian perturbation theory (which is common practice in 
grid based simulations). This baryon density field is initially 
Gaussian. 
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ILPT/LLA: Again ILPT is used, but now the initial baryon 
density field is computed according to equation (35), so that 



it is consistent with ILPT. This initial baryon density field is 
no longer Gaussian. 

2LPT/LLA: Second order Lagrangian perturbation theory 
(cf. Section 3.1.2) is used for dark matter and baryons, us- 
ing equation (35) a consistent initial baryon density field is 
generated. 

We find that a similar transient behaviour as that between 



ILPT and 2LPT for pure CDM simulations (cf. eg. Crocce 
et al.||2006| [Tatekawa Mizuno 2007) can be observed. The 
skewness of the density field is underestimated when starting 
at Zs — 100 with the initially Gaussian field (EPT). This is 
somewhat improved when starting at higher redshifts. When 
using an initial density field that is consistent with ILPT, the 
skewness is boosted so that the dependence of the PDF on the 
initial redshift of the simulation becomes smaller. When us- 
ing 2LPT and LLA, no significant dependence of the PDF on 
the starting redshift can be observed. Note that we show the 
results for the baryon density field directly taken from the sim- 
ulations rather than a density field that was smoothed with an 
aperture filter. Filtering with a finite apperture however leads 
to almost identical results, the skewness is underestimated at 
high redshifts independent of scale. 

In summary, Eulerian perturbation theory leads to an un- 
derestimation of the skewness of the baryon density PDF at 
high redshifts. When using the LLA for the baryons, a PDF is 
imposed that is consistent with Lagrangian perturbation the- 
ory. In particular, as has been shown already for dark matter 
simulations, the density PDF becomes much less sensitive to 
the initial starting redshift when using 2LPT. Using the local 
Lagrangian approximation, consistent initial baryon density 
fields can be construed. 



6 SUMMARY & CONCLUSIONS 

We have presented and implemented a method to generate 
real-space sampled initial conditions for nested grids, extend- 
ing prior work by [Pen ( 1997 ) , [Bertschinger (2001) and Sirko 
('2005). These initial conditions are suitable for "zoom-in" cos- 
mological simulations of structure formation with several lev- 
els of refinement that allow the study of single cosmological 
objects, such as e.g. groups and clusters of galaxies, single 
galaxies, or first stars, in a cosmological context at very high 
resolution while maintaining a small computer memory foot- 
print. 

Apart from the general advantage of real-space sampled 
initial conditions discussed by Pen ( 1997) and "Sirko (2005), 



namely a correct reproduction of real-space statistical proper- 
ties, such as the two-point correlation function and mass vari- 
ance in spheres, we demonstrate that our multi-scale convolu- 
tion approach for nested grids has very favourable error prop- 
erties: (1) no interference of different modes between different 
refinement levels, and (2) errors confined to coarse-fine bound- 
aries. We find RMS errors of velocity/ displacement fields in 
the refined region at the order of 10~^ in units of the stan- 
dard deviation of the respective fields which is an improvement 
of about two orders of magnitude over previous approaches. 

In order to determine particle displacement and velocity 
fields, at first and second order Lagrangian perturbation the- 
ory, we followed two approaches. First, a pure multi-grid Pois- 
son solver is used. We show that the finite volume discretiza- 



tion it assumes leads to a suppression of perturbations at the 
smallest scales that can however be corrected by a deconvo- 
lution. Second, a hybrid Poisson solver is developed, which 
uses an adaptive multi-grid algorithm for inter-level gravity 
and an FFT-based Poisson solver for the finest grid. This hy- 
brid approach leads to no suppression of small scale perturba- 
tions. Analyzing statistical properties of unigrid simulations 
for which initial conditions have been generated with the two 
different approaches, we find however that non-linear objects 
are not sensitive to the lack of power at these very small scales 
and differences are only seen in the mildly non-linear regime 
and at the lowest halo masses. We conclude that our real- 
space based approach is well suited also for initial conditions 
for unigrid simulations. 

In order to verify the accuracy of nested initial conditions 
generated with our method, we investigated the properties of 
a galaxy cluster both in a unigrid simulation and a "zoom- 
in" simulation with one and two levels of refinement. We find 
that the gross properties of the cluster, such as virial mass, 
radius, spin, shape and velocity dispersion, are recovered at 
per cent level or better in the re-simulations. Density profiles 
show no bias with scatter at the level of a few per cent, mainly 
attributable to slight changes in the positions of sub-structures 
in the re-simulations. We also find that the sub- halo mass 
function is recovered with an accuracy of a few per cent in 
all re-simulations, and observe that the hybrid Poisson solver 
performs slightly better than the pure multi-grid approach. 
We thus conclude that our algorithm does not introduce bias 
or unreasonable scatter in observables that will be deduced 
from such a re-simulation and will thus provide a reliable tool 
to study the internal structure of cosmological observables at 
high resolution in large cosmological volumes. 

Finally, we study the inclusion of a baryonic component 
in our approach by generalising the use of real space transfer 
functions also for distinct baryon and CDM density and veloc- 
ity perturbations. We demonstrate that our approach repro- 
duces the evolution expected from linear perturbation theory 
correctly. Furthermore, we propose to set the initial baryon 
density field based on a local Lagrangian approximation which 
is consistent with Lagrangian perturbation theory of first or 
second order and which greatly reduces the dependence of the 
baryon density field on the starting redshift of the simulation, 
thus reducing transient behaviour. 
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